Label-free liquid biopsy through the identification of tumor cells by machine learning-powered tomographic phase imaging flow cytometry

Image-based identification of circulating tumor cells in microfluidic cytometry condition is one of the most challenging perspectives in the Liquid Biopsy scenario. Here we show a machine learning-powered tomographic phase imaging flow cytometry system capable to provide high-throughput 3D phase-contrast tomograms of each single cell. In fact, we show that discrimination of tumor cells against white blood cells is potentially achievable with the aid of artificial intelligence in a label-free flow-cyto-tomography method. We propose a hierarchical machine learning decision-maker, working on a set of features calculated from the 3D tomograms of the cells’ refractive index. We prove that 3D morphological features are adequately distinctive to identify tumor cells versus the white blood cell background in the first stage and, moreover, in recognizing the tumor type at the second decision step. Proof-of-concept experiments are shown, in which two different tumor cell lines, namely neuroblastoma cancer cells and ovarian cancer cells, are used against monocytes. The reported results allow claiming the identification of tumor cells with a success rate higher than 97% and with an accuracy over 97% in discriminating between the two cancer cell types, thus opening in a near future the route to a new Liquid Biopsy tool for detecting and classifying circulating tumor cells in blood by stain-free method.

www.nature.com/scientificreports/ above decision, with the aim to collect only the identified CTCs. The idea follows a conceptual workflow similar to the so-called Intelligent Image-Activated Cell Sorter (iIACS) 49 . The results attained here show a tumor cells identification success > 97 % along with an accuracy > 97 % in discriminating between the two cancer cell types. Such challenging capability confirms, in principle, the feasibility of the proposed method in realistic scenario of future LB by stain-free mode.

Results
Tomographic dataset and 3D data augmentation. The TPI-FC experimental setup and the corresponding numerical processing described in the "Methods" section have been used for the experiments and for retrieving the reported tomographic dataset. The dataset has been used for training a hierarchical ML model at the aim of identifying first tumor cells from MCs and then recognizing the type of tumor. In particular, hundreds of digital holograms of each single cell were recorded while the cells were flowing and rotating along a microfluidic channel. From each digital hologram, the corresponding quantitative phase map (QPM) is numerically extracted. In a QPM, the information about the 3D spatial distribution of the cell's RIs is coupled to the information about the 3D morphology of the cell in the form of a 2D image. Therefore, to decouple them, after retrieving the corresponding unknown rolling angles, the QPMs are combined to reconstruct the 3D RI tomogram at the single-cell level. Following this pipeline, the 3D RI tomograms of several cell lines have been reconstructed, as summarized in Table S1. In particular, (1) as regards the MC class, 247 THP1 cells have been Steps of the machine learning-powered tomographic phase imaging flow cytometry. (i) A holographic video sequence containing cells while flowing and rotating along a microfluidic channel is recorded; (ii) for each cell, hundreds of QPMs are numerically computed from the recorded digital holograms; (iii) for each cell, the several QPMs are used to reconstruct its 3D RI tomogram; (iv) the 3D RI tomogram is given in input to a hierarchical ML classifier to detect whether it is a monocyte or a tumor cell and, in the second case, whether it is a neuroblastoma cancer cell or an ovarian cancer cell. www.nature.com/scientificreports/ recorded, (2) as regards the NB class, 372 cells have been recorded, i.e. 115 CHP212 cells, 106 SKNBE2 cells, and 151 SKNSH cells, (3) as regards the OC class, 310 cells have been recorded, i.e. 95 A2780 cells and 215 CAOV3 cells. In Fig. 2a,c,e, an example of QPM is reported for a MC cell (THP1), a NB cell (SKNSH), and an OC cell (CAOV3), respectively, while the central slices of the corresponding 3D RI tomograms are shown in Fig. 2b,d,f, respectively. For these three cells, the overall sequence of QPMs and their 3D RI reconstructions (both a sliceby-slice visualization and an isolevels representation) can be seen in the Supplementary Movie S1 . To train and  test a ML model, the tomographic dataset has been divided into a training set containing 700 cells (200 MC,  250 NB, and 250 OC cells) and a test set containing 229 cells (47 MC, 122 NB, and 60 OC cells), as reported in  Table S1. To counteract the limited sizes of our datasets, the training set has undergone a 3D data augmentation process. In particular, 9×, 3× and 3× augmented tomograms have been computed from each MC, NB, and OC cells respectively. Therefore, as shown in Table S1, the augmented tomographic dataset is made of 2000 MC, 1000 NB, and 1000 OC cells. To perform the 3D data augmentation, three successive operations have been applied to each 3D RI tomogram, i.e. (1) intensity scaling, (2) intensity shifting, (3) morphological alteration. The 3D data augmentation, as well as the holographic processing, the tomographic reconstructions, the feature extraction, and the ML training have been carried out by means of MATLAB® R2022b. Let n x, y, z be the RI volumetric distribution inside the L x × L y × L z array containing the cell (in our TPI-FC system, L x = L y = L z = 201 pixels), and let n 0 be the RI of the surrounding medium, supposed homogeneous. Let Ŵ be the volumetric support of the cell, i.e. the set of voxels x, y, z occupied by the cell, and let n Ŵ x, y, z be the set of cell's RI values. In our experiments, n 0 = 1.334 and �n Ŵ x, y, z > 0 , where �n Ŵ x, y, z = n Ŵ x, y, z − n 0 .  www.nature.com/scientificreports/ 3. As regards the morphological alteration, the L x × L y × L z volume is rescaled into the (L x + c x ) × L y + c y × (L z + c z ) volume through the imresize MATLAB function, where c x , c y and c z are separately randomly drawn from the uniform distribution U(−20, 20) , thus obtaining the �n x, y, z values ( Ŵ is the new volumetric support of the cell after the morphological alteration). Hence, for example, if c x > 0 , the tomogram is stretched along the x-axis, otherwise, if c x < 0 , the tomogram is compressed along the x-axis.
After these three operations, the augmented 3D RI tomogram is obtained as n x, y, z = �n (3) x, y, z + n 0 . An example of intensity scaling, intensity shifting, and morphological alteration is reported in Fig. 2g-i, respectively, based on the 3D RI tomogram of the SKNSH cell in Fig. 2d.
Features extraction. In order to train a ML model, 44 features have been calculated from each 3D RI tomogram. In particular, 11 features are related to the RI statistical distribution: the average, the median, the mode, the maximum, the standard deviation, the skewness, the entropy, the kurtosis, the 0.25-quantile, the 0.75-quantile, and the dry mass (defined as the mass of the cell without considering its aqueous content 50 ). Other 9 features are instead related to the 3D morphology of the cell, such as the volume, the convex volume, the sphericity, the extent, the solidity, the lengths of the three principal axes, and the normalized centroids' distance. The convex volume is the volume of the smallest convex polygon containing the cell. The sphericity is the ratio between the surface area of a sphere with the same volume of the cell and the surface area of the cell (it is 1 for spherical objects, otherwise it is less than 1). The extent is the ratio between the cells volume and the volume of the cells bounding box, i.e. the smallest cuboid containing the cell. The solidity is the ratio between the volume and the convex volume. The principal axes are the major axes of the ellipsoid having the same normalized second central moments as the cell. The normalized centroids' distance is the distance between the centroid and the weighted centroid (with respect to the RI volumetric distribution) of the cell, normalized to the cells equivalent radius, defined as the radius of a sphere having the same volume as the cell. The morphological features have been computed by means of the regionprops3 MATLAB function. Finally, 24 features have been computed as the Haralick features of the 3D Gray-Level Co-Occurrence Matrix (GLCM) 51 . The GLCM describes the different combinations of the grey levels within an image. Indeed, the GLCM G(h, k, i x , i y , i z , d) measures how many times a pixel with value h occurs along the direction i x , i y , i z at distance d with respect to a pixel with value k . Therefore, after fixing 12 Haralick features have been computed 52 , i.e. energy, entropy, correlation, contrast, variance, sum average, inertia, cluster shade, cluster tendency, homogeneity, maximum probability, and inverse variance. The 12 Haralick features at d = 0.5 μm have been finally averaged among the 13 3D directions, thus obtaining 12 of the 24 GLCM features. Instead, the other 12 GLCM features have been computed in the same way by fixing d = 1 μm. The histograms of some features computed from the 3D RI tomograms of the augmented training set are displayed in Fig. 3a-l. In particular, in Fig. 3a-f, the average RI, the volume, the dry mass, the sphericity, the GLCM energy, and the GLCM contrast are shown by separating the MCs from the cancer cells (NB and OC cells). Instead, the same features are reported in Fig. 3g-l, respectively, by separating the NB cells from the OC cells. For example, it can be inferred that the sphericity histogram of the tumor cells is left-shifted with respect to that of the MCs (see Fig. 3d), while the average RI histogram of the OC cells is right-shifted with respect to that of the NB cells (see Fig. 3g). Furthermore, in the following Section, results of the ML classification based on the 3D RI tomograms will be compared with the ML classification based on the 2D QPMs of the same cells. Therefore, 44 features like those described above have been measured from each QPM. The 11 features related to the phase statistical distribution are again the average, the median, the mode, the maximum, the standard deviation, the skewness, the entropy, the kurtosis, the 0.25-quantile, the 0.75-quantile, and the dry mass 28 . The 9 2D morphological features are the area, the extent, the solidity, the circularity, the eccentricity, the maximum and minimum Feret diameters, the major axis length, and the normalized centroids' distance, computed by means of the regionprops MAT-LAB function. The circularity is computed as 4Aπ/P 2 , where A and P are the area and the perimeter of the cell, respectively, and it is 1 for a perfect circle. The major axis length is the length of the major axis of the ellipse having the same normalized second central moments as the cell, while the eccentricity is the ratio of the distance between the foci of the same ellipse and the major axis length. The maximum and minimum Feret diameters are respectively the maximum and minimum distances between any two boundary points on the antipodal vertices of the convex hull enclosing the cell. The normalized centroids' distance is the distance between the centroid and the weighted centroid (with respect to the phase spatial distribution) of the cell, normalized to the cell's equivalent radius, defined as the radius of a circle having the same area as the cell. Hierarchical ML classification. Usually, in LB, CTCs are firstly detected inside a background of WBCs, and then the particular type of cancer is recognized. Therefore, following this scheme, we have built a hierarchical classifier made of two successive binary classifiers (namely, A and B). A shallow neural network with three layers has been chosen for both classifiers A and B, as sketched in Fig. 4a. The input layer is made of 44 nodes, corresponding to the 44 features described in the previous Section. The hidden layer is a fully convolutional layer made of 100 nodes for the classifier A and 10 nodes for the classifier B. The output layer is made of 1 node  4c and Table 1). However, within the hierarchical classifier, the misclassification error propagates from the input of the classifier A to the output of the classifier B. Hence, after extracting the 44 features from the 3D RI tomograms of the single-cells recorded through the TPI-FC system, they are given as input of the classifier A, which detects a MC or a tumor cell with a certain error. The 44 features of the sole cells classified as tumor cells at the output of the classifier A are given as input of the classifier B, which at the end recognizes the type of cancer (i.e. NB or OC) with its additional error. Therefore, the hierarchical classifier detects the MC cells with REC = 95.7% , the NB cells with REC = 93.4% , and the OC cells with REC = 98.3% (Fig. 4d and Table 2). To measure the overall performance of a classifier, we considered the accuracy (ACC), that is the probability of correctly classifying an element given as input of the model. The classifier A has ACC = 97.4% (Fig. 4b and Table 1), the classifier B has ACC = 97.3% (Fig. 4c and Table 1), and the overall hierarchical classifier has ACC = 95.2% (Fig. 4d and Table 2). Hence, the combination proposed here between the TPI-FC data and a ML hierarchical classifier offers the possibility for recognizing and then phenotyping cancer cells with very high accuracy. This is mainly due to the ability of TPI-FC of providing a statistical measurement of the distinctive fingerprint of the single cells belonging to a certain population (thanks to the FC module) because their biophysical content can be accessed in 3D (thanks to the TPI module). However, the tomographic reconstruction requires larger hardware resources and longer computational times if compared to a more conventional QPI in flow cytometry (QPI-FC) system, because hundreds of digital holograms per cell must be processed instead of just one digital hologram per cell. For this reason, a comparison of the classification performance achievable by a QPI-FC system is needed to justify the greater burden of a www.nature.com/scientificreports/ TPI-FC system. To this aim, to have a fair comparison between these two strategies, by using the straight-ray model of the light propagation (i.e., the same model used by the employed tomographic reconstruction code), one QPM per cell has been computed by integrating at 0 • each 3D RI tomogram belonging to the augmented training set and to the test set described in Table S1. Then, the hierarchical classifier has been trained by means of the 44 features measured from the 2D QPMs. In this case, the logistic regression has been chosen as classifier A and the linear discriminant has been chosen as classifier B, since they provided the best test performance (see Tables 1 and 2). In particular, the classifier A has ACC = 87.3% (instead of the 3D ACC = 97.4% , see Table 1), the classifier B has ACC = 85.2% (instead of the 3D ACC = 97.3% , see Table 1), and the overall hierarchical classifier has ACC = 76.0% (instead of the 3D ACC = 95.2% , see Table 2). Moreover, as reported in Table 2  www.nature.com/scientificreports/ REC = 98.3% ). In summary, all the performance given by a QPI-FC system are drastically lower than the TPI-FC method, as summarized by the accuracy of the hierarchical classifier, which decreases by 19.2%.

Discussion
The proposed methodology based on the ML-powered TPI-FC has been demonstrated to be very efficient in identifying tumor cells with respect to the background of MCs, also providing the possibility to discriminate between two cancer cell lines, i.e. NB and OC. While only a proof-of-concept study has been addressed in this paper, it is reasonable to foresee that the extension of the approach to the identification of other cancer cell lines, as well as other types of WBCs, is feasible. In this case, it is needed to consider a more heterogeneous population of WBCs and CTCs. Moreover, for the cancer cells, further levels of classification could be included for the phenotyping of CTCs, thus discriminating among sub-types. On the other hand, the choice of using only MCs as the WBC population has two main motivations. First, it has been demonstrated that MCs may play an important role in metastasis, and they can extravasate and differentiate into macrophages, promoting tumor cell extravasation, survival, and subsequent growth 53 . Therefore, the possibility to analyze MCs in combination with cancer cells can provide further information about tumor treatments. Second, MCs are often comparable in size with respect to CTCs, thus making a challenge their preliminary separation from the whole blood sample. In fact, since the intent of the proposed approach would be the integration in existing technologies that provide the removal of platelets, red blood cells and smaller WBCs, there is high probability to collect MCs and CTCs together. Our results show that ML-powered TPI-FC can correctly identify tumor cells with respect to MCs with an accuracy > 97 %, thus extremely limiting this scenario. In addition, the possibility to use the classification output as the actuator signal in an intelligent image-activate sorter further confirms the ability to increase the CTCs enrichment effectiveness. However, to clearly understand the potentiality of the proposed method in a real LB application, it would be valuable to quantify the device performance using well-recognized parameters 54 such as capture efficiency, enrichment, purity, throughput, cell viability, and release efficiency. By considering the experiments carried out in this paper, we can define the purity as the probability of correctly recognizing a tumor cell with respect to a background of interfering cells, i.e. the MCs, that is 97.8 % (i.e., the recall parameter). The cell viability is expected to be very high in our label-free system, and cells can be used again for further www.nature.com/scientificreports/ downstream analyses, mainly represented by high-throughput sequencing approaches, that allow us searching in unbiased and large-scale manner for different cancer genomic alterations, or approaches for mutation specific detection. As regards the throughput, we have already demonstrated that, for a Field-of-View (FoV) = 160 µm × 215 µm, the holographic data for reconstructing an upper bound of about 50/60 tomograms/s can be recorded 41 .
Hence, in the current FoV = 640 µm × 640 µm, this value can be increased of about 10 times. Moreover, the release efficiency, which describes the number of cells that are recovered by the system, should be related to the ability of the system in providing a video-rate processing and classification of the tomograms. Finally, the capture efficiency and enrichment, which are related to the ability of the system in capturing CTCs from a sample, cannot be defined at this stage since, as stated above, we have considered already done the pre-filtration of the blood sample. In summary, identification of CTCs, reflecting intra-tumoral heterogeneity, represent a promising pathway for cancer diagnosis and monitoring although the current unavailability of universal and specific cell-surface markers makes CTCs detection a challenging task. In this context, we have proposed a ML-powered TPI-FC system aimed to discriminate tumour cells within a blood cells background and to catch their phenotypic heterogeneity. Our results suggest that this label-free approach for rapid and efficient phenotyping of cancer cells could represent a suitable tool for the identification of novel morphological biomarkers to discriminate CTCs and to distinguish clinically aggressive from more favourable CTCs. TPI-FC system. The opto-fluidic recording system sketched in Fig. 5a has been employed for the TPI-FC experiments 40 . The optical module is an off-axis DH microscope based on a Mach-Zehnder interferometer. The beam generated by a coherent laser source (Laser Quantum Torus emitting at 532 nm) is divided into an object and a reference beam by a Polarizing Beam Splitter (PBS). Two Half-Wave Plates (HWPs) are used to balance the ratio between the two beams' intensities. After illuminating the cells flowing inside the microfluidic channel, the object beam is collected by a high numerical aperture microscope objective (MO1, Zeiss Plan-ApoChromat 40×, NA = 1.3) and sent to a tube lens (TL1). Then, a Beam Splitter (BS) takes in input the object beam and the reference beam, which has passed in the meantime through a beam expander, a second microscope objective (MO2) and a second tube lens (TL2). The BS produces the interference between the two beams, properly collimated by the TLs, thus creating the digital hologram recorded at 30 fps by the CMOS camera (Teledyne Dalsa Genie Nano-CXP 5120 × 5120 pixels, 4.5 µm/pixel, 80 fps) and stored by digital recording unit (IOindustries DVR, Core 2, Video Storage Modules 4 × 3840 GB). The fluidic module is composed of an automatic low-pressure pump (CETONI base 120 Nemesys) generating a laminar flow at about 50 nl/s inside a rectangular-cross-section microfluidic channel (Microfluidic ChipShop 10,000,107 -200 μm × 1000 μm × 58.5 mm). Therefore, cells flow along the microchannel and contemporarily rotate due to the torque generated by the velocity gradient in the channel cross-section, which, in turn, is due to the parabolic velocity profile characterizing the laminar flow. Let z be the optical axis and y be the flow direction. Cells are mainly recorded in the center with respect to the x -axis and in the bottom with respect to the z-axis 52 . As the FoV measures 640 µm × 640 µm (5120 × 5120 pixels), hundreds of digital holograms per cell are recorded at 30 fps along different viewing angles. An example of full-FoV digital hologram is displayed in Fig. 5b. Despite the excellent performance achieved here, the throughput of the device can be further increased by optimizing the design and the flow conditions of the microfluidic module, which is the object of ongoing work. www.nature.com/scientificreports/ Holo-tomographic numerical processing. Each 5120 × 5120 digital hologram of the video sequence recorded by the TPI-FC system is numerically processed to obtain the QPMs of flowing single cells 28 . As shown in Fig. 5b,c, from each holographic frame, 384 × 384 region of interests (ROIs) are cropped around the cells. As www.nature.com/scientificreports/ the DH microscope is in off-axis configuration, each holographic ROI is demodulated by selecting the real diffraction order in the Fourier spectrum through a band-pass filter 56 . In DH, the focal plane of each cell can be recovered numerically after the experiment. As displayed in Fig. 5c, the holographic ROI (i.e., the 384 × 384 red box) is larger than the size of the single cell (i.e., the 201 × 201 dashed green box) in order to include its entire diffraction pattern. It is precisely the diffraction pattern that allows implementing an automatic strategy for refocusing the cell (namely, autofocusing) 57 . In fact, the cell's diffraction pattern changes when the demodulated holographic ROI is numerically propagated along the optical z-axis through the Angular Spectrum method 56 . For each z-distance, this variation can be quantified as the image contrast of the complex wavefront's amplitude, measured through the Tamura Coefficient (TC). The focal distance of each cell is computed by minimizing its TC functional 57 . Then, the argument of the refocused complex wavefront is extracted, residual optical aberrations are removed by subtracting a reference hologram without cells 58 , and the resulting phase map is denoised through the 2D windowed Fourier transform filtering 59 and unwrapped through the PUMA algorithm 60 . At the end of this process, a 201 × 201 QPM is selected around the cell's weighted centroid 37 , as reported in Fig. 5c,d. Therefore, each 384 × 384 holographic ROI is converted into the corresponding 201 × 201 QPM. To reconstruct the 3D RI tomogram of a single cell, the hundreds of QPMs recorded during its roto-translation along the microfluidic channel are combined together through the Filtered Back Projection (FBP) algorithm 29 , which takes in input also the corresponding viewing angles. In the TPI-FC system, the viewing angles correspond to the rolling angles of the cell, which are estimated by exploiting the microfluidic properties and the holographic tracking 55 , since cells mainly rotate around the x-axis with a continuous and quasi-uniform roto-translational speed.

Methods
In a conventional fluorescence imaging flow cytometry system, usually just one 2D intensity map per cell is recorded. Instead, in our TPI-FC system, hundreds of digital holograms per cell are acquired, which are used to reconstruct its 3D RI tomogram after computing the corresponding 2D QPMs. By using an interpreted programming language (i.e., Matlab ® 2022b) over an Intel ® Core™ i9-9900K CPU with a 64Gb RAM, the numerical processing to compute a QPM from the holographic ROI takes 7.71 s per frame. This numerical processing must be replicated for all the hundreds of recorded holographic ROIs related to each single cell, thus representing at this stage the actual bottleneck for reconstructing very large datasets of 3D RI tomograms (in fact, the final tomographic reconstruction based on the FBP algorithm takes just few seconds per cell). However, the holographic processing can be speeded up by 45 times through deep learning, thus reducing the computational time of the phase retrieval from 7.71 s to 0.17 s per frame, but at the cost of losing the finest details inside the QPMs 40 . Moreover, another way to reduce the computational burden is to consider that reconstruction of all the QPMs of one single cell, used to obtain its 3D RI tomogram, is a highly parallelizable process. Thus, multicore and/ or multimachine implementations of the reconstruction software could run in parallel on GPUs. Instead, the computational time for detecting the holographic ROIs belonging to one single cell is much lower with respect to the computational time for their QPMs retrieval. In particular, the detection of all the cells inside each frame of the holographic sequence (see Fig. 5b) takes 0.16 s. The following process for grouping all the holographic ROIs belonging to the same cell takes about 20 s per cell. For example, for a cell running the entire FoV in 100 frames, the first detection step takes 100 × 0.16 s = 16 s, which must be added to the 20 s of the second detection step. However, it is important to specify that, while the second step is performed separately for each cell, the first step is shared by all the cells flowing along the channel at the same moment, as shown in Fig. 5b, thus it is implemented just once for all of them.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.